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ABSTRACT 

Hydrodynamic models are becoming prevalent design tools for small scale devices and other 
devices in which high energy effects can dominate transport. Most current hydrodynamic models use 
a parabolic band approximation to obtain fairly simple conservation equations. Interest in accounting 
for band structure effects in hydrodynamic device simulation has begun to grow since parabolic 
models can not fully describe the transport in state of the art devices due to the distribution populating 
non-parabolic states within the band. This paper presents two different non-parabolic formulations of 
the hydrodynamic model suitable for the simulation of inhomogeneous semiconductor devices. The 
first formulation uses the Kane dispersion relationship <hk)*/2m - W(1 + aW) . The second 
formulation makes use of a power law {(Wt)V2m - xW*} for the dispersion relation. Hydrodynamic 
models which use the first formulation rely on the binomial expansion to obtain moment equations 
with closed form coefficients. This limits the energy range over which the model is valid. The power 
law formulation readily produces closed form coefficients similar to those obtained using the parabolic 
band approximation. However, the fitting parameters (x,y) are only valid over a limited energy range. 
The physical significance of the band non-parabolicity is discussed as well as the 
advantages\disadvantages and approximations of the two non-parabolic models. A companion paper 
describes device simulations based on the three dispersion relationships; parabolic, Kane dispersion, 
and power law dispersion. 
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Carrier Energy 

Planck's constant divided by 2n 
Reciprocal lattice vector 

Non*parabolicity factor under the Kane dispersion relation, positional dependent 

Carrier mass at the band edge, a constant but positional dependent 

Adjustable parameters for the power law dispersion relation, both positional dependent 

Density of states in momentum space 

Conduction band edge 

Electron quasi-fermi level 

Boltzmann's constant 

Electron temperature 

Fermi integral of various order 

Fermi integral divided by Gamma function 

Gamma function 

Effective density of states 

Reduced energy {(E fc -E„)/(KT,)} 

Semiconductor band gap 
Gradient operator in physical space 
Gradient operator in k space 
Group velocity 
Distribution function 
Equilibrium distribution function 
Relaxation time 

Conduction band potentials (electrostatic, affinity, and band gap narrowing) 

Electron concentration 
Identity matrix 
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INTRODUCTION 


The use of hydrodynamic models for device simulation are becoming common as 
characteristic device dimensions continue to decrease. Hot electron effects can play a dominant role 
in carrier transport for high mobility semiconductors. Current hydrodynamic models consist of a set 
of conservation equations derived by taking moments of die Boltzmann transport equation (BTE). 
During the derivation of the conservation equations the parabolic band approximation is used to obtain 
rather simple coefficients on the forcing terms in die flux equations. By relying on die parabolic band 
approximation higher order energy transport effects due to variation in the band structure are 
neglected. Interest in accounting for band structure effects in hydrodynamic device simulation has 
begun to grow because parabolic models can not adequately account for high energy effects in 
semiconductors with non-parabolic band structures. Non-parabolic band formulations have a history 
dating back to the 1950's [1-3]. However drift-diffusion models and more specifically hydrodynamic 
simulators with non-parabolic band formulations are a very recent topic of research. 

Several non-parabolic hydrodynamic models have been reported for homogeneous material 
systems [4-7] using the Kane dispersion relationship [3]. The general functional form obtained is 
similar to parabolic hydrodynamic models with first order corrections on the diffusion term. Azoff [8] 
derived a hydrodynamic model suitable for degenerate heterostructure semiconductors though the final 
form of the equations was not directly amenable to current device simulation codes. However, Azoff 
clearly showed that a forcing term due to a gradient in the non-parabolicity factor exists. Woolard et. 
•1- [9] presented a non-parabolic hydrodynamic model based on moments of the velocity and energy 
(n, W(k)) instead of the momentum and momentum squared (hk, tiV). This leads to a simpler energy 

conservation equation. However, the non-parabolic coefficient in the field term and the forcing terms 

* 

due to non-uniform band structure were neglected in the other moment equations. Cassi and Riccb 
[10] introduced an alternative to the Kane relation in the form of a power law for the dispersion 
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relationship. Instead of using the classical Kane dispersion law relating the energy and momentum, 
die band was fit over a specified energy range using two adjustable parameters. The approximations 
and assumptions implied by assuming the power law formulation were absent. It will be shown below 
dial the power law formulation of the dispersion relation leads to a more simplistic and compact 
formulation than the classical Kane expression. 

The purpose of this paper is to introduce two non-parabolic hydrodynamic formulations 
suitable for the simulation of devices with inhomogeneous material layers. The final form of the 
conservation equations will be in a form which will allow incorporation into existing device simulation 
codes, similar to the parabolic formulation. These different formulations are based on different 
choices for the dispersion relationship, one uses the standard Kane dispersion for non-parabolic bands 
and the second uses a power law relationship [10]. For comparison the parabolic hydrodynamic 
formulation will also be presented. The form of the conservation equations are strongly affected by 
the non-parabolicity factor of the bands, the choice of the dispersion relationship, and the assumptions 
made to simplify the coefficients. As in the case of the parabolic formulation, both non-parabolic 
formulations require estimates of higher order moments to provide mathematical closure of the 
relationships. It will be shown that more physical insight can be obtained by examining the terms 
from the power law formulation due to their similarity to the parabolic formulation. The power law is 
advantageous in that the terms of the conservation equations are the same as in the parabolic 
formulation, except for a simple multiplicative constant when Boltzmann statistics are employed. 


Dispersion Relations and Carrier Concentration 
The two non-parabolic dispersion relations relating the energy to the momentum are 



Where a is the non-parabolicity factor in the Kane dispersion relation. This parameter is usually 
calculated from a k*P perturbation approach to the band structure. The a factor can also be estimated 
in terms of known parameters of the semiconductor [10]. The second non-parabolic dispersion relation is die 
power law formulation of Cassi and Ricci [10] which is not a first order approach to the band structure as 
compared to the Kane dispersion relation. Instead, the parameters x and y are obtained by a best fit to 
equation (la) over a specified energy range. In reference [10] this energy range was (1.3 eV, 3.0 eV) and 
the authors produced a very good fit for a - 0.4789, determined by inference from x and y and is reproduced 
in Figure 1. The formulae for the earner concentration using these two dispersion relations and assuming 
Fermi-Dirac statistics are (the binomial expansion has been employed for equation (2a)) 
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In the case of the Kane dispersion the non-parabolicity is a perturbation to the parabolic model 
regardless of the fact that Fermi-Dirac statistics were used in die derivation instead of Boltzmann 
statistics. In the power law case, unlike the a formulation, die carrier concentration will only differ 
from the parabolic by a multiplicative constant if Boltzmann statistics are assumed. 

Within the hydrodynamic simulation of a device the carrier concentration is calculated by 
directly solving the continuity equations. Inspection of equations (2a) and (2b) shows that the 
concentration in both cases is a function of the reduced energy, X\. The determination of 11 at each 
position is crucial to the successful simulation of a device since most of the factors within the 
transport equations i.e. number of ionized dopants, Shockley-Read-Hall recombination etc., are 
functions of T] as well. Therefore, from the knowledge of the carrier concentration at any specific 
position the reduced energy, T|, is calculated by suitably inverting equations (2a) and (2b). With these 
facts in mind it is illustrative to compare the non-parabolic formulations for the earner concentration 
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to flic parabolic formulation at various reduced eneigy values. Figure 2 displays the deviation from 
the parabolic formulation for a - 0.4789, x - 1 365, and y - 1.52 [10]. The case of the frill Kane 
dispersion relation before binomial expansion is also given in Figure 2 and shows that the binomial 
expansion is justified for this case of the non-parabolicity parameter, the two curves are almost 
identical throughout the entire reduced energy range. Notice also that the deviation of the alpha 
nonparabolicity formulations is at most 5% from the parabolic case until the reduced energy exceeds 
zero, when the difference in the orders of the Fermi integral is larger. However, in the case of the 
power law fit at high energy, the deviation from the parabolic case is more than 82% and more than 
87% from the a formulation. To explain this large difference the plots of the dispersion relations must 
be re-examined. Figure 1. At the high energy part of the curve the power law and a formulations 
have nearly identical dispersion relations. The insert in Figure 1 shows the low energy range and 
shows the maximum deviation of the two curves. Though there is a greater deviation between the two 
curves at lower energy than at high energy, the maximum difference between the two curves is 
insufficient by itself to fully account for the large difference in the carrier concentrations. However, 
carrier concentration depends upon both the density of states and the distribution function. Since the 
distribution function weights more heavily to lower energy states, a mismatch at low energy is 
strongly reflected in the product of the distribution function times the density of states and 
subsequently the carrier concentration. Therefore, in order to obtain a matched value of the carrier 
concentration, it is critical to have a closer fit to the low eneigy range of the band. For comparison, 
a low energy fit was made using a different set of for x and y, x - 1.185 and y - 1.052. These values 
of x and y were determined by fitting the power law dispersion relation over the low eneigy range 
(0.0 eV^W^ 0.2 eV). As can be seen from the inset of Figure 1, the dispersion relations match 
very closely at low energy but deviate significantly at higher energy. Figure 2 also shows the 
calculated carrier concentrations for various reduced energy values using the low energy power law fit 
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factors. Compared to the power law fit to the high energy range, using the power law fit to the low 
energy range of the band produces concentration values which closely match those obtained using the 
a formulation. 

Particle Flux Equations 

The general flux (0 - v,) conservation equation is determined by taking moments of the BTE 

as 
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Notice that the third, fourth and fifth terms on the left hand side are tensor products. The factor of 3 
in die fifth term is due to the order of parenthesis in the original moment equation. This moment 
equation can not be processed further until some functional form of the effective mass is assumed, 
which depends on the choice of the dispersion relation. Using the parabolic dispersion and die two 
dispersion relations in equation (1), the resulting flux equations are (positional subscript on the 
gradient operators has been removed) 
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Each of these equations must be integrated over all k apace or equivalently over energy wing the 
density of states [15]. Before the integration is performed one more assumption must be made, drat 
the relaxation time is independent of k or W. If the constant relaxation time assumption is not made 
then the energy dependence of the relaxation time must be moved through the gradient operator on the 
second term in each of equations (4-6) and a term accounting for the gradient of the relaxation time 
must be re-created. Changing the integration from k space to energy space, substituting for the 
mobility, Gi-r/m), and making v/ 0 equal to zero gives the following integral equations 
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To produce closed form solutions for the integrals b equation (8) the bbomial expansion is used 
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repeatedly and all terms of order a? or higher are aet to zero, to finally become 
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The only term in the power law flux equation (equation (9)) which will require an expansion is the Vy 
term due to the ln(W) factor in the integrand. The expansion may only be required for certain choices 
of the distribution, but to maintain generality it is applied for all distributions. A parabolic 
interpolation, using the points 0, 0.5, and 1, to the part of the equation containing the log term is 
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With the above substitution the power law flux equation can be written as 
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The approximation int equation (1 1) tends to degrade as y increases. However, this forcing term will 
not appear in simulations of homogeneous materials. In the case of inhomogeneous materials systems 

this term may be quite small as y has only a limited range, (1 £ y £ 2). 

One of the goals of this paper is to derive a hydrodynamic model suitable for fast efficient 
simulation of state of the art devices. As previously stated the parabolic model is inadequate for 
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certain material systems due to the fact that their band structures deviate dramatically from a parabolic 
shape at high energies. In these situations, a non-parabolic band structure provides a more accurate 
description. As discussed above, the full a formulation captures die desired physics, but is 
unattractive because the coefficients within die transport equations need to be numerically evaluated. 

In contrast, the binomial a formulation and the power law formulation provide closed form 
expressions for the transport coefficients. However, the question remains as to their ability to display 
the proper trends; how closely do the coefficients match the full a formulation? We will show that 
the binomial a formulation has a limited energy and non-parabolicity range due to the binomial 
approximation. Non-physical results will be obtained if the formulation is extended into regions 
outside the binomial limit. The power law does not produce non-physical results but more closely 
matches the parabolic formulation when fit to the low energy part of the band. 

The coefficients on similar forcing terms in the various formulations can now be compared to 
examine the impact of the non-parabolicity factor and approximations made during the derivation. 

The comparison is done at this time to avoid any confusion from assuming a distribution function. 
Table I lists the forcing terms which are compared and the terms within the integrands involved. 

Table II contains the values of the non-parabolicity factors at which these factors are compared, this 
includes values of a from 0.04 to 4.0 and appropriate (xy) values fit to these O values over two 
different energy ranges (0.0, 0.2) and (1.5, 3.0). Figure 3 shows the prefactor for the diffusion term 
from all the flux equations; Figure 3a is for slightly non-parabolic bands (a - 0.04); 3b is for the 
non-parabolicity in reference [10] (a « 0.4789), and 3c is for a highly non-parabolic band (o - 4.0). 
Figures 4 a< and 5 a-c display the prefactors on the forcing terms from changes in the effective mass 

and the field for the three degrees of non-parabolicity. For all three forcing terms at the lowest values 
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of non-parabolicity (Figures 3a, 4a, 5a) the prefactors compare favorably with the parabolic 
formulation. For the diffusion term (Figures 3 a-c) as the non-parabolicity factor increases the two 
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cases of the power law formulation match very closely to the full a formulation. On die other hand 
the a formulation that utilizes the binomial expansion is very different, especially at the highest value 
of non-parabolicity considered. From Figure 3c it is clear that die binomial a formulation is clearly 
incorrect if the energy exceeds 0.3 eV. Even before this point the diffusion will be underestimated. 

In die case of the mass term. Figure 4, die power law formulation which is fit over the low energy 
range is much closer to the parabolic case, as expected due to die small change in the fitting 
parameters from their parabolic values. However, the power law with the parameters fit over a larger 
energy range more closely matches the full a formulation. As in die case of the diffusion term, the 
binomial a formulation severely underestimates the effect of this forcing term especially as the non- 
parabolicity is increased, and is limited to energies less than 0.1 eV for a ■ 4.0. The coefficients for 
the field term. Figure 3, follow the same conclusions as for the mass term. From these figures it is 
clear that the binomial a formulation has a very limited energy range of validity as the non- 
parabolicity factor is increased. Using this formulation at higher energies or high non-parabolicity 
factors can give un-physical results due to the prefactors changing sign. On the other hand, the power 
law formulation with parameters fit over a small energy range will tend to produce results which more 
closely match the parabolic band model. It will not produce un-physical trends and does appear to 
have a larger range of validity for both energy and non-parabolicity factors. The case of the power 
law with parameters fit over a large energy range more closely matches the full a formulation in terms 
of the forcing coefficients. However, due to the problems previously described for the calculation of 
the carrier concentration this advantage may be immaterial. 

Table I can also be used to gain some physical insights into transport in the non-parabolic 
band structures, especially using the power law formulation. First notice that as the non-parabolicity 

t 

factor is decreased (a-+0, x— »1, y-tl) all the coefficients reduce to die parabolic case indicating that 
all three formulations are equivalent in this respect. As the non-parabolicity factor and energy 
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increases the binomial a formulation can actually predict a change in the sign of a forcing term. In 
the full a formulation the sign on the forcing terms does not change but it is unclear as to how the 
term decreases. In the case of the power law, when the parameter y equals 2 it is obvious dial the 
field term will become identically aero. The only terms which will be non-zero in die flux equation 
will be the diffusion, gradient in mass, and gradient in the non-parabolicity terms. This can be 
explained with the use of the power law energy equation (equation (lb)} and the group velocity 
equation. When y equals 2 there is a linear relationship between die momentum and energy, die 
bands are V shaped. The group velocity is proportional to the gradient of the energy with respect to 
k, which for the case of y equals 2 means that the group velocity is a constant. Therefore, no matter 
how much force is applied to the election by external forces the velocity is not increased. The only 
factors which can produce current are gradients which can change this fixed velocity; i.e. changes in 
mass or non-parabolicity factor, or changes in the number of earners moving at this fixed velocity, 
diffusion. Therefore, non-parabolic formulations which do not include coefficients to diminish the 
field term as the non-parabolicity increases overestimate the flux. 

At this point in the derivation a recursion relation must be formulated, a distribution function 
assumed, or some other mathematical method (Minimum-Maximum theorem) must be used to provide 
mathematical closure for equations (7), (10), and (12). Since, a goal of this paper is the formulation 
of models suitable for the numerical simulation of devices the first two options are explored (16). A 
recursion relation would allow moments of higher order to be approximated by lower order moments, 
the lower orders are calculated from the conservation equations. This option does not require that a 
specific form of the distribution function be used, unless the recursion relations are based on a specific 

distribution. However, in the case of non-parabolic bands the standard recursion relations may no 

* 

longer be applicable (9). Therefore, this option was not pursued. The other option, and the one 
chosen for this work, is to assume a specific form for the distribution function, higher moments can 
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then be calculated based on the known distribution function. Some of die choices for the distribution 
function include heated Maxwellian, shifted and heated Maxwellian, heated Fenni-Dirac, or shifted 
and heated Fenni-Dirac. Since the Maxwellian distributions can be recovered by relaxing die 
degeneracy, the Fenni-Dirac distributions were the only ones considered for this work. 

In the non-parabolic formulation a simple relation between energy and velocity will not exist 
due to die change in die density of states. Also higher order powers of the energy are required to 
close the relationships in the a formulation, this will require cross product terms involving die 
temperature and the velocity. In addition, the power law formulation has non-integer powers of the 
energy which will be very difficult to evaluate for the shifted and heated distributions. Due to these 
conditions and the fact that all the formulations break down as the energy rises, the heated Fenni- 
Dirac distribution was used to close the relationships. The flux equation in the binomial a formulation 
is 
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The flux equation in the power law formulation becomes 
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The reader can verify that as the non-parabolicity factors are diminished, (a-X>, the two non- 

parabolic formulations reduce to the standard parabolic case. The equations also reduce to simpler 
forms when the degeneracy effects are ignored (all orders of Fermi integral reduce to exponentials). 
The flux equations can be discretized using normal techniques to produce comparable simulation 
codes under various assumptions [17]. 

Energy Flux Equations 

For the energy flux equation (© - Wv,) the general conservation equation is 
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Notice that the fourth through seventh terms on the left hand side are tensor products. As in the case 
of the particle flux moment the energy flux moment equation can not be processed further until some 
functional form of the effective mass is assumed. By making similar assumptions, substitutions, and 
approximations {binomial expansion, Fermi -Dirac statistics, equipartition of energy ...} die energy flux 


equations for the three dispersion relationships become 
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First, equation (16) corrects a sign error on the order of a Fermi integral which occurred in an earlier 
publication [17]. Again, the reader can verify that the non-parabolic formulations reduce to the 
parabolic case as the non-parabolicity factors are decreased. As in the case of the particle flux 
equations the energy flux equations can be discretized using normal techniques [17]. 


CONCLUSIONS 

Two formulations of the hydrodynamic model have been presented for the simulation of non- 
parabolic inhomogeneous material systems, the standard Kane formulation and the power law 
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formulation of Cassi and Riccd. Both forms reduce to the parabolic hydrodynamic model as the non- 
parabolicity factors are diminished. The forcing terms in the particle and energy flux equations due to 
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variations in the non-parabolicity factors of the hands as suggested by AzofiT have been recreated in 
both non-parabolic formulations. It was shown that the binomial a formulation is suitable for the 
calculation of the carrier concentration but has a limited energy and non-parabolicity range when 
applied to die coefficients of the flux equations. Extending the binomial a formulation past these 
limits leads to non-physical terms in both the particle and energy flux equations. In the case of the 
power law formulation it was shown that when the adjustable parameters were fit to the high energy 
range the deviation in the carrier concentration from both the parabolic and O cases was extreme. 

When fit to the lower energy range the power law produced carrier concentrations comparable to the a 
formulations. The energy and non-parabolicity range of the power law formulation for the particle 
and energy flux equations is larger than the binomial a formulation, but it is still limited by die 
adjustable parameters. However, unlike the binomial a formulation die power law will not lead to 
physically unrealistic results, but will tend to more closely match the parabolic formulation when the 
adjustable parameters are fit over the low energy portion of the band. It is shown and argued using 
the power law formulation and a high non-parabolicity factor (y - 2) that a dampening factor must 
exist on the field term in the particle flux equation to account for the non-parabolicity of the bands in 
both non-parabolic formulations. Particle flux and energy flux equations using a heated Fermi-Dirac 
distribution function are also presented, to allow for incorporation into existing device simulators. 
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Table I. Terms (prefactors) within the integrands of the forcing coefficients which are graphically 
compared. 


Formulation 

Equation 

Number 

Diffusion Term 

Mass Term Vm 

Field Term Ve„ 
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Table II. Values of the non-parabolicity factors which are compared. This includes O values from 
0.04 to 4.0 and (x,y) values fit to these a values over two different energy ranges (0.0, 0.2) and (1.5, 

3.0). 


Case 

a 

x,y {0.0, 0.2} 

x,y {1.5, 3.0} 

a 

1 

1.0148, 1.0045 

1.0207, 1.0827 

b 

0.4789 

1.185,1.052 

1.363, 1.52 

c 

4.0 

2.975,1.322 

4.8233, 1.901 
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Figure Captions 

Figure 1 Dispersion relations for the Kane model (solid line, a - 0.4789), the power law formulation 
(dashed iine, x-1.365, y - 1.52 fit over high energy, 1.5 S W * 3.0 eV) and the power law 
formulation (dotted line, x-1.185, y - 1.052 fit over low energy, 0.0 £ W £ 0.2 eV). The insert 
shows the detail of the low energy range. 

Figure 2 Deviation, from die parabolic case, of the carrier concentration as a function of the reduced 
energy for the full O formulation (solid line, O - 0.4789), binomial a (dotted), power law fit at high 
energy (dashed, x«l 365, y-1.52), and power law fit at low energy (dashed, x-1.185, y-1.052) 
dispersion relations as a function of the reduced energy (t)) at T,— 300 K. 

Figure 3. Comparison of the integrand coefficients (prefactors) occurring on the diffusion term of the 
particle flux formulations. See Table I and H for the exact form of the equations and degrees of non- 
parabolicity considered. 

Figure 4. Comparison of the integrand coefficients (prefactors) occurring on the gradient of the 
effective mass term of the particle flux formulations. See Table I and II for the exact form of the 
equations and degrees of non*parabolicity considered. 

Figure 5. Comparison of the integrand coefficients (prefactors) occurring on the gradient of the 
potential term occurring in the particle flux formulations. See Table I and II for the exact form of 
the equations and degrees of non-parabolicity considered. 
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